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We study the evolution of linear density fluctuations of free-streaming massive neutrinos at redshift 
of 2 < 1000, with an explicit justification on the use of a fluid approximation. We solve the 
collisionless Boltzmann equation in an Einstein de-Sitter (EdS) universe, truncating the Boltzmann 
hierarchy at Zmax = 1 and 2, and compare the resulting density contrast of neutrinos, 5^""^, with 
that of the exact solutions of the Boltzmann equation that we derive in this paper. Roughly 
speaking, the fluid approximation is accurate if neutrinos were already non-relativistic when the 
neutrino density fluctuation of a given wavenumber entered the horizon. We find that the fluid 
approximation is accurate at few to 25% for massive neutrinos with 0.05 < m,y < 0.5 eV at the 
scale of fc < 0.4 h Mpc~^ and redshift of 2 < 10. This result quantifies the limitation of the fluid 
approximation, for the massive neutrinos with ^0.5 eV. We also find that the density contrast 
calculated from fluid equations (i.e., continuity and Euler equations) becomes a better approximation 
at a lower redshift, and the accuracy can be further improved by including an anisotropic stress term 
in the Euler equation. The anisotropic stress term effectively increases the pressure term by a factor 
of 9/5. 

PACS numbers: 98.65.Dx 98.70.Vc 98.80.Cq 



I. INTRODUCTION 

What is the mass of neutrinos? We know that at least 
two of three standard model neutrino species have finite 
masses. The constraints on the squared mass differences 
of the three species of neutrinos obtained from solar 
fisl and atmospheric oscillation experiments [T^-fTs} (re- 
views can be found in [l9l42^ ) are (3 cr errors) 



eV^ 



Ami, - (7.65t°:^^o) x 10"^ 
|A<| = (2.40t°i^)xlO-3 eV^ 



(1) 
(2) 



Therefore, the lower limit on the sum of neutrino masses 
is 0.058 eV. Observations of the CMB and large-scale 
structure of the universe can provide limits on the ab- 
solute mass of neutrinos. The current upper bounds on 
the sum of neutrino masses are ~ 0.3 — 0.6 eV j23l - l28j . 
In this paper, we use 0.58 eV (95% CL) from WMAP7yr 
as a conservative upper bound on the sum of neutrino 
masses [29l |. 

The large-scale structure of the universe is a sensitive 
probe of neutrino masses [l^, [30l - l39| . Massive neutrinos 
suppress the small-scale matter power spectrum by their 
large velocity dispersion. The fractional amount of the 
suppression in the small-scale limit is roughly given by 



AP(fc) 



with 



P(fc) 



94.1 eV 



(3) 



(4) 
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where the summation is taken over the i-th species of 
neutrinos [H, EO, SJ . 

Relativistic neutrinos are not a fiuid. Massive neu- 
trinos, being collisionless, are not a fiuid, either. How- 
ever, when the velocity dispersion of massive neutrinos 
becomes low enough, they may be approximately treated 
as a fluid, just as we normally treat Cold Dark Matter 
(CDM) particles as a fluid on large-scales. While this 
is a reasonable expectation, as far as we know, the ex- 
tent to which the fluid approximation is valid for massive 
neutrinos has not been discussed in the literature. 

Then, why is a fluid approximation useful while we 
have Boltzmann codes such as CMBfast [42| and CAMB 
[43j . which solve the Boltzmann equations numerically 
to the accuracy of order ^ 0.1%? First, we will have 
more physical insight to the growth of the neutrino den- 
sity fluctuations by directly solving continuity and Euler 
equations rather than numerically solving a set of infinite 
order of Boltzmann hierarchy. Second, and most impor- 
tantly, as the density fiuctuations become non- linear, i.e., 
(5 1, we need to use higher-order perturbation theories 
to accurately model the small-scale density fiuctuations. 
Since the higher-order perturbation theories have been 
constructed for CDM with a fluid approximation, we can- 
not simply modify theories to include massive neutrinos if 
a fluid approximation is not valid for those particles. On 
the other hand, if a fluid approximation is valid for some 
range of redshifts, length scales and neutrino masses, we 
can greatly simplify the model of non-linear density fluc- 
tuations in the presence of massive neutrino, as shown in 

In this paper, we shall study the validity of a fluid ap- 
proximation of massive neutrinos. To achieve this goal, 
we first solve the Boltzmann equations describing the 
evolution of the perturbed phase-space distribution func- 
tion of massive neutrinos exactly and compare the exact 



2 



results to the results with the fluid approximation, i.e., 
solutions with the higher multipole moments (/ > 3) ig- 
nored. Then, we shall examine the ranges of applicability 
of fluid approximation in both spatial and time scales, as 
a function of neutrino masses. 

The rest of this paper is organized as follows. In § 
we briefly review the effects of massive neutrino free- 
streaming on the structure formation of the universe. In 
§ mil we provide the basic fluid equations and the lin- 
earized Boltzmann equation required for our theoretical 
flame work. In § IIVI we briefly discuss the analytic so- 
lutions of the Boltzmann equation for collision-less par- 
ticles. In § |Vl we compare the exact solutions of the 
Boltzmann equations with the fluid approximation, and 
discuss the limitation of the fluid approximation for sev- 
eral masses of massive neutrino. Finally, in § I VII we 
discuss the implications of our results and conclude. In 
Appendix]^ we discuss how to define the free-streaming 
scale starting from the fluid equations. In Appendix |BJ 
we give the detailed derivation of the exact solution of 
the Boltzmann equation both for massless and massive 
neutrinos. Even though our main interest is in massive 
neutrinos, our results shown here are also applicable to 
collision-less particles in general, whose time evolution of 
the perturbed phase space distribution follows the lin- 
earized collision-less Boltzmann equation with the zero- 
th order phase space distribution function being frozen at 
sufficiently early time (i.e., we set the initial conditions 
of the neutrino transfer function after the decoupling of 
neutrino, ~ 1 MeV). 



II. THE FREE-STREAMING OF THE MASSIVE 
NEUTRINO 

We are interested in the mass range of 0.05 < m^^j < 
0.58 eV for the most massive species of neutrinos, which 
became non-relativistic well after the matter radiation 
equality. The mass density of the massive neutrinos rel- 
ative to the total matter density is given by 



1 



n^h^ n,nh^ 94.1cV ' 



(5) 



where the summation is taken over the different species 
of neutrinos. Neutrinos become non-relativistic when the 
mean energy per particle of neutrinos in the relativistic 
limit, 



(E) 



J d^pp {cxp[p/T,{z)] + l)-^ 
Jd^p (exp[p/r,(z)] + l)-i 
7tt^ 



180C(3) 



(6) 



falls below m^j. By solving 3.15r^,o(l + ^nr) = m-jy^i, one 
finds the redshift of relativistic to non-relativistic transi- 
tion epoch, Ziir, as 



(7) 



for the i-th neutrino species. 

The density fluctuation of neutrinos cannot grow 
within the horizon size until neutrinos become non- 
relativistic. Once neutrinos become non-relativistic, the 
neutrino density fluctuation begins to grow on scale 
greater than the so called "free-streaming scale," which 
is set by the velocity dispersion of neutrinos: 



Jd^p pVmg i(exp[p/r.(z)] + l) 

Jd^p (exp[p/T,(z)] + 
15C(5) f4_\ ^ Tl,{l + 2)2 
C(3) [u 



(8) 



where p is the pro per momentum of the massive neutrino 
(see Appendix of |4l|). 

The wavenumber corresponding to the free-streaming 
scale, fcpSj is defined by the single- fluid continuity and 
Euler equations: 



<5(k,T) +0(k,r) ^0 

e{k,T) + niT)9{k,T) + 



(9) 

(5(k,T) = 0, 
(10) 



where 



0.677 



~ V 2 cs(z) 



3 n{z) 

2 o-i/,i(z) 



(Sv) i^^n{i + z)' + n, 



h Mpc"^ 
(11) 



Here, derivatives are with respect to a conformal time, 
dr — dt/a, %{t) = and 6'(k, r) is a velocity di- 

vergence of the fluid. Note that Eq.(l8|) assumes that 
neutrinos are non-relativistic. 

In Figure [U we show /cps.i from Eq.(fTT|) (dotted line), 
comoving horizon scale, aH{a), (thick solid line) and 
/cps.i calculated numerically from Eq.®, where m^^j is 

replaced by ^^p^ + j (thin solid line). In this figure, 

we use m,y.i = 0.13 cV. 

We find that the free-streaming scale is close to the 
horizon size until the relativistic to non-relativistic tran- 
sition of a neutrino, and once the neutrino becomes 
non-relativistic, the free-streaming scale decreases as 
A:fs(<^) ^ Let ^3 examine the evolution of the neu- 

trino density fluctuations at three length scales: 



^ Here, we say Cs ~ cr^ i; however, strictly speaking, the velocity 
dispersion defined in Eq.||8ll should not be used to define the free- 
streaming scale, fcps , as the Euler equation contains sound speed, 
Cg = not the velocity dispersion. In the non-relativistic 

limit, we have Cs = ^g^o"^,i — 0.745cr^_i. We derive this relation 
in Appendix 1X1 
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FIG. 1: Free-streaming scale of a massive neutrino, feps.i, (black line), comoving horizon scale, aH{a), (thick black line) and an 
approximation to the free-streaming scale in the non-relativistic limit given by Ea. (|lip . (dotted line) as functions of the scale 
factor, a. We use mi^j = 0.13 eV. The horizontal lines show (1) large, (2) small, and (3) intermediate scale modes as described 
in §11] 



1. At the large-scale, where k ^ kps{a) for all a < qq 
(ao is the present-day scale factor), the neutrino 
density fluctuation starts to grow soon after the 
mode enters the horizon, and its time evolution is 
identical to that of CDM, 6,^{k,a) = 6c{k,a). 

2. At the small-scale, where k ^ kpsio-) for all a < uq, 
the neutrino density fluctuation oscillates around 
its initial value due to the free-streaming effect, 
St,{k, a) ~ 6„{k, Oi) ~ 0. 

3. At the intermediate scale, the mode first experi- 
ences the free-streaming phase, and thus does not 
grow. Once k < kpsia) is satisfied, the mode starts 
to grow, rapidly catching up with the gravitational 
potential set up by CDM. 



III. THE BOLTZMANN HIERARCHY AND 
FLUID APPROXIMATION 

In this section, we provide all the relevant equations 
and definitions needed for our theoretical flame work, 
following [iBj in the conformal-Newtonian gauge. 

For fermions and bosons, we have the phase space dis- 
tribution (in natural units) given by 



/o(g,r) 



5s 



(27r)3 ee('?,r)/aT(a) ^ 



(12) 



The linear order perturbation to the distribution func- 
tion, vl'(k, n, q, r), is defined as 

/(k,n,Q,T) = /o(g,r)[l + \E'(k,n,(z,T)], (13) 

where q = |q| and n = q/q. 

Since neutrinos decoupled while they were highly rel- 
ativistic, the unperturbed distribution function after the 
neutrino decoupling continues to be given by its relativis- 
tic form: 



/o(g) 



1 



(27r)3 e9/'^'^(") ± 1 ' 



(14) 



even after neutrinos become non-relativistic. The tem- 
perature of such collision-less particles decreases as 
T(a) = To{ao/a), even when they are non-relativistic. 

The evolution of the linearized phase-space distribu- 
tion for collision-less particles such as CDM and neutri- 
nos is governed by the linearized collision-less Boltzmann 
equation. 



d^{k, n,q,T) 



dr 



+ i 



(k ■ n)*(k, n,q, t) 



dln/o(g) 
dlnq 



0(fc,r) -i^^^(k-n)^/.(fc,T) 



= 0, (15) 



where ip and (j) are a Newtonian gravitational potential 
and a curvature perturbation, respectively. ^ 



where the sign of "+" is for fermions and " is for 
bosons, q and e(g, r) = \/ q'^ + a'^{T)m? are the co- 
moving momentum (i.e., q = a(T)p) and the comoving 
energy of a particle, respectively. Here, r is a con- 
formal time, which is related to the proper time by 
dr = dt/a{t), and gs is a number of degrees of freedom. 



In the original work of [45ll . ip and <f> are defined as scalar 
perturbations in the metric in the conformal Newtonian gauge: 
ds2 ^ „2(^j[_(;^ _|_ 2ip)dT^ -I- (1 - 24,)dx* dxi]. They are related 
to the gauge invariant variables and of |4a | and 'I' and 
3> of [43] hy tp = <3>A = 'I' and </> = — <I>£f = — <I>. 
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To simplify the equation, we define \f(k, n, g,T) s 

\l/(k, n, r) (^^^^1^^^ , and replace the time deriva- 
tive from T to X = fcr, and re- write Ea. (|15p as 



a*(k,n,g,x) . q ~ , d<t>{k,x) 

+ t-, rAi* n, q, x) + 

ox e[q, X) ax 

-i^ill^^^(^k,x) ^0, (16) 



where ^ is a cosine between the wavenumber and momen- 
tum, i.e., k • n = kfi. Finally, we expand the Boltzmann 
equation (Eq. (jl6p ) by Legendre polynomials, using 



*(k,n,g,x) = 5^(-z)'(2; + l)*,(fc,<z,a;)Pz(Ai),(17) 



and obtain a set of infinite series of differential equations 
(also known as Boltzmann hierarchy) as follows: 



%{k,q,x) 
^[{k,q,x) 

^'l{k,q,x) 



e{q,x 

q 



3e((?, a;) 
3q 



[^oik,q,x) ~2^'2ik,q,x)] 

ij{k,x), (19) 
[;*/_i(fc,g,x) 



(2Z + l)e(g,a;)' 
- il + l)^i+iik,q,x)] (for;>2), (20) 



where the primes denote derivatives with respect to 
X = kr. Here, 'iio{k,q,x) is sourced by ^i{k,q,x). All 
the successive multipoles with / > 1, '^i>i{k,q,x), are 
sourced by ^i-i{k,q,x) and 5';+i(fc, g, x), so that the 
evolution of Z-th multipole propagates the whole system 
of equations back and forth. In order to close the sys- 
tem of equations, we need to truncate the Boltzmann 
hierarchy at some finite multipole, Zmax- Now, the ques- 
tion is, "m which condition the fluid approximation (i.e., 
^max ~ 1 or 2) is valid?^ 

To make a contact with the familiar form of fluid equa- 
tions, we relate multipoles of the perturbed distribution 
function, 5';(fc, g, r), to the quantities such as the density 

contrast, (5(fc,r) = ^^y^, velocity dispersion, 0{k,T), 
and anisotropic stress, a{k,T), by integrating '^i{k,q^T) 



over the momentum space with appropriate powers of q 
in 



P{r) 
P{r) 
dp{k,T) 
6P{k,T) 
{p+P)0ik,T) 
{p+P)a{k,T) 



An 
3^4(7) 

47r 



q dq e((7,T)/o(g), 

^2 



q dq 



■fo{q), 



(21) 
(22) 



q'dq eiq,T)fo{q)^o{k,q,T), (23) 



lq^dq-^fo{q)Mk,q,r)m 



3a4(r) 
Ank 



q^dqqUq)^,ik,q,T), (25) 



^'^ I q^dq-^Mq)^2{k,q,r). {26) 



3a4(r) 



We obtain the fluid equations by truncating the Boltz- 
mann hierarchy at Imax = 2: 



<j(fc,T) = -[l + w{T)][0{k,T) -30(A:,r)] 



SP{k,T) 



w(r) 



e{k,T) 



Sp{k, T 

[l-3w{T)]0{k,T 



(27) 



_a{T) 
a{r) 



5ik,T), 

w{t) 
1 + w{t) 



e{k,T) 



1 -I- w{t) 



(28) 



■^l{k,q,x)~<t)'{k,x), (18) ^(fc,T) 



°^^^-[2-3u;(T)](T(fc,T) !i^ilI^a(fc,T) 



a(r)' 

-e(fc,r) + 4iE(fc,r), 
15 a{T) 



1 + w(t) 



(29) 



where w{t) = P{t)/p{t) is an equation of state, and we 
have defined the following variables: 



(p + P)e(fc,r) 



(p + P)S(fc,T) 



Ank I 2, I q 
q dq q' 



a\T)J ^ ^ \eiq,T) 
X /o(g)*i(fc,g,r). 



(30) 



87r 



9 rf?- 



3a4(r)y ""e(g,T) V^l^.r), 
X fo{q)^2ik,q,T). (31) 



In the relativistic and non-rclativistic limits, where ma- 
jority of neutrinos in the phase space distribution have 
momenta of g ^ ^{q-,''') and 9 ^ ^{q^T): '^e have 
(w,w, f|,e,E) = (i,0, i, 61, cr) and (0,0,0,0,0), respec- 
tively. 

Since CDM is non-relativistic throughout the redshift 
of our interest, we can greatly simplify the calculation of 
the density contrast of CDM by fluid approximation (i.e., 
Zmax = 1). We have 



S{k,T) = -e'(fc,r) +3</.(fc,r). 



0{k,T) 



a(r) 
a(r) 



6'(fc,T) +fcV(fc,T). 



(32) 
(33) 
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As for massive neutrinos, the fluid approximation may 
or may not be valid, depending on the mass, scale or 
redshift of interest. We check the validity of the fluid 
approximation for massive neutrinos by comparing to the 
exact solutions in section § |Vl 

When I max — 2, Eg. (|18p and (|20|) give a useful relation 
between '^o{k,q,x) and '^2{k,q,x): 



d_ 

dx 



2 ~ 2 
*2(fc, g, x) + -*o(A:, g, x) + -(t>{k, x) 
5 5 



0, (34) 



which gives 



2 ~ 2 
^2{k, q, x) + --^aik, q, x) + -(f>{k, x) 

2 ~ 2 
= fzCfc, q, Xi) + -*o(fc, g, x^) + Xj), (35) 

where a;^ ^ 1 is an initial time. With this relation and 
Eas. (p4| and p6|. we can rewrite the anisotropic stress, 
a{k,T), in the Euler equation (Eg. ([28])) in terms of pres- 
sure, SP{k,T), as 



k^a{k,T) + - 



4SPik,T)/dp{k,T) 



1 + w{t) 



k'^d{k, t) = const, 



(36) 



where we have set (j) = const. At late times, t ^ Ti, 
where S{k,T) :s> S{k,Ti) and (j{k,T) ^ a{k,Ti), the right 
hand side of Eg.([36|) is negligible compared to the second 
term on the left hand side. Therefore, we have 



fc2cr(/fc,r) 



ASPik,T)/Sp{k,T),^ 



1 + w(t) 



k^5{k,T) 



(37) 



This result shows that a increases the pressure by a factor 
of |. The free-streaming scale, fcps: is then reduced by a 
factor of ~ 1.34. 



IV. ANALYTIC SOLUTIONS FOR THE 
BOLTZMANN EQUATION 

In this section, we briefly describe analytic solutions of 
the Boltzmann eguation (Eg. ([16])), to which the fluid ap- 
proximation is compared. We give a detailed derivation 
of the solutions in Appendix |B] 

Instead of expanding the Boltzmann eguation by Leg- 
endre polynomials as in [45|, we first find a formal solu- 
tion of Eg.([l| 



^ik,q,ii,x) = ^ik,q,fi,x,)e-'^^'<^^^-'<^^'^^ 

+ [ (ix'e-'^l"(^)-^(^')lS'(A:,g,//,x')(38) 



where 



Qii ^ - -£(9'^) in \ 9(p{k,x) 

S{k,q,p,x) = 1 p'4>{k,x) , 39 

Q ax 



and wc have defined 

Z{x) E 



e{q,x') 



-dx'. 



(40) 



The lower integration boundary, C G 7?., is an arbitrary 
constant. We then expand the above formal solution with 
Legcndrc polynomials. We will need the following rela- 
tion of the Wigner 3-j symbols, 

^P,(^)P,(^)P„(^,) = (^J J . (41) 

The solution for ^'/(A:, q, x) for a given I is given by 
Mk, q,x)^Yl E(-*)''^'""'(2/' + l)(2r' + 1) 

2 



(' I" 



x-i>i'{k,q,Xi)ji"{z - Zj) 







ip{k) / dx' 



,e{q,x') 



I 



21 + 1 



ji-i{z-z') 



1 



21 + 1 



(42) 



where wc have assumed ip{k,x) = 0(fc,x) = (which 
is satisfied in an Einstein de-Sitter, EdS, universe). We 
derive solutions for the most general case (i.e., ip{k, x) ^ 

and (pikjX) ^ 0) in Appendix IB] 

As we see in Eg. ip^ . the infinite scries of the Boltz- 
mann hierarchy, Eqs. ([T51) ~ (pni) . is now expressed in 
terms of the spherical Bessel functions, ji{x), its inte- 
grals weighted by e{q,x)/q, and the infinite sum of the 
initial values of 5';(fc, q, x). 

While Ea. (|^^ appears to have infinite sums over /, the 
sum actually truncates because, at initial time (which is 
taken to be before the horizon re-entry), '^i{k,q,x) for 

1 > 3 can be ignored H^. Together with the triangular 
inegualities of the Wigner 3-j symbols. 



\l-l'\ < I" <l + l' 



(43) 



only finite terms remain in the solution of '^i{k,q,x). 
The explicit solutions of ^'o(fc, q, x) and 4'i(fc, q, x) are 



*o(fc, 9, x) ^o{k, q, Xi)jo{z - Zi) 
-3^i(fc, g,.T,;)ji(z - Zi) + 5'^2ik,q,Xi)j2{z - Z.i) 
,t{q,x') 



+%l){k) I dx' 



-ii{z - z'), 



(44) 



^i{k,q,x) = ^ii{k,q,Xi)ii{z - Zi) 
+'^i{k,q,Xi)jo{z ~ Zi) " 2*2(^,9, Xi)ji(z - Zj) 
-2*i(fc, g,.T,;)j2(2 - Zj) +3*2(^,(7,2:^)^3(2; - Zj) 



Mk) 



dx' 



,e{q,x') 



jo{z - z') - -j2{z - z') 



(45) 



Let us examine the behavior of 5'i(fc, q, a;). At suf- 
ficiently late time, z(x) ^ z(xi), all the terms con- 
taining initial values of 4';(fc,q, x) become negligible, as 
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ji{z) — > for z ^ The last term, which does not 
depend on initial values, is the dominant term. For rela- 
tivistic neutrinos, e{q, x) — q, the last term is propor- 
tional to dx'ji(x — x'), which approaches constant 
for x ^ I. For non-relativistic neutrinos, e{q,x) = 
a{T)m cx x'^, ji{z) does not oscillate for a; ^ 1 (because 
-^(^) = / ^(^'^2; cx i). Thus, the integrand, ^^ji(z), 
grows. 

V. THE VALIDITY OF THE FLUID 
APPROXIMATION 



we need to calculate the '^[){k^q,x) for relativistic par- 
ticles including all the higher multipoles in principle. If 
the neutrinos arc sufficiently massive, then the last term 
in Eq. (|19p becomes dominant over the first two terms 
for most of the particles in the phase-space distribution, 
fo{q)i making the fluid approximation (i.e., Zmax = 1 or 
2) valid. 

Note that a fluid approximation does not imply that 
all the higher moments of the Boltzmann equations are 
small. It just means that the evolution equations of 
S{k,x) and 6{k,x) are decoupled from the higher mo- 
ments of the Boltzmann equations. 



Before we start, let us remember why fluid approxima- 
tion may be valid for massive particles. Ea. (|19|) shows 
that the I = 2 mode becomes unimportant when the 
gravitational force (the last term) becomes dominant. 
The ratio of the last term to the first two terms is of 
order (e/q)^, which is unity for relativistic particles, but 
it is much greater than unity for non-relativistic particles. 
Thus, I > 2 modes become irrelevant for the evolution of 
I = and 1 modes, allowing us to truncate the Boltz- 
mann hierarchy at Zmax = 1- 

For example, when the distribution function of neutri- 
nos is dominated by the non-relativistic states (i.e., q <^ 
a(x)m), we have the following solutions for 4'o(fc, x) 
and 5'i(fc,q, a;) with constant and ip^ 



^o(K,g,a;) = — 



(46) 

where C is a constant, and the fastest growing 
modes grow as \^o(fc,(7, x) cx a(x) and \E'i(fc, 5, a;) cx 
a{x) a/1 -I- a{x). 

The observable quantities such as i5(A:,a;) and 9{k,x) 
are given by the integrals of ^o{k,q,x) and ^i{k,q,x). 
They pick up contributions from the relativistic particles 
{e{q,x) ~ q) as well, but the phase-space number density 
of those relativistic particles is exponentially suppressed. 
To see this, we rewrite /o of relativistic particles as 



/o(g) 



(2^) 



3 



± 1 



(48) 



where ^ > 1 for relativistic particles in the phase space 
distribution. 

For CDM with m ^ 1 GeV, the mass to temperature 
ratio, m/aT{a), is very large for the time that is relevant 
to the structure formation; thus, only extremely non- 
relativistic particles, q/m <^ 1, contribute to the den- 
sity contrast, S{k, r). As neutrinos are much lighter than 
CDM, relativistic particles may or may not contribute to 
the density contrast significantly. For example, to cal- 
culate the massless neutrino density contrast accurately. 



A. ^i{k,q,x) with various Zmax 

To check the validity of fluid approximation, we need to 
compare the exact solutions of ^0(^1 ^) a-nd 4'i(fc, g, x) 
to the approximate solutions of 4'o(fc, 9, x) and 4'i(fc, g, x) 
with finite ^max = 1:2 and 3. To simplify the problem, 
we make the following (reasonable) approximations: 

1. The universe is flat and matter dominated (EdS), 
for which <f> = ip = Q. 

2. We ignore the evolution of and caused by 
massive neutrinos, which is a good approximation 
for proving the validity of the fluid approximation. 
That is, V' is determined by CDM only. 

3. The initial perturbations are adiabatic and the 
wave lengths of the initial perturbations are greater 
than the horizon size, i.e., Xi <^ 1. Speciflcally, we 
have (Eq.(97) of [H): 



^o{k,q,Xi) = -^(5i.(fc,Xj) = 



'^i{k,q,Xi) 



1 e[q,Xi) 9^{k,Xi) 
3 q k 

1 e{q,x^) 



6 q 



-Xitp, 



(49) 



(50) 



^2{k,q,Xi) = -'^cr^ik,Xi) ^ -^xftjj, (51) 



*,>3(fc,g,.T,) = 0, 



(52) 



where the most of neutrinos in the phase-space distribu- 
tion are initially highly relativistic, e(g,.T,;) ^ q. 

We rewrite the Boltzmann equations fEas. (fT8| ~ (jSO]) ) 
in terms of dimensionless parameters k/C and m/q de- 
fined below: the ratio of the comoving energy to the co- 
moving momentum of a particle is given as 



^ = Wl + a^(x)(^) , 



and the scale factor is given as 

, , C fCx 



(53) 



(54) 
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where 

c = y = H,^n^y—j 

= 0.0183Vf^ h Mpc"\ (55) 

with Coq = 1, Peg = p(acq), and we assume the mat- 
ter radiation equality to happen at 1 + Zoq = ao /ocq = 
3000. For Vtr = 8.47 x IQ-^^ a» = 0-25 and C = 
0.0092/1 Mpc-^ 

With this convention for the scale factor, given comov- 
ing momentum is equal to the physical momentum (i.e., 
q = p) ai the matter radiation equality, a = aeq, and 
therefore, Ea. (l48)) becomes, ^ 

(27r)'^ g-i^irr ± 1 



We see that the asymptotic growth rate is 
^oxactj-^^ ^ 2;) oc oc a (see Eq.fgn])). For neutri- 
nos with m/q = 10 and k/C = 100, the growth of 
yj;cxactj-^^ ^-j -g suppressed between the horizon re- 
entry and the epoch of relativistic to non-relativistic 
transition [xnr = 8.3); however, once neutrinos become 
non-relativistic, '^'^^'^'^{k,q,x) grows rapidly, catching 
up with the gravitational potential set up by CDM. ^ 
For both large and small-scales, fluid approximation 
becomes more accurate as we increase the /max, but for 
neutrinos with m/q — 100 and k/C = 10, /max = 1 is 
sufficient to approximate the exact solution to better 
than 1% accuracy for almost entire evolution history 
{xh < X < xq). For small-scale, k/C = 100, with nearly 
relativistic neutrinos {m/q = 10), fluid approximation 
with low multipole cutoff (i.e., /max = 1,2 and 3) breaks 
down for almost entire evolution history, while at late 
time (a; > 1000), we start to see an accuracy of better 
than 1% for /max = 2 and 3. 



1. ■^o{k,q,x) 

Figure [2] shows the evolution of ^o{k,q,x) for two 
different scales {k/C ~ 10 and 100) and two differ- 
ent momenta {m/q = 10 and 100). We calculate 
^q"'^(A:, g, x) by truncating the Boltzmann equations at 
^max = 1,2 and 3 (fluid approximation), while we calcu- 
late \I'o™'^*(/c, g, x) from the exact solution of the Boltz- 
mann equations given in Ea. (|44p . The fractional er- 
ror in the fluid approximation is defined as A^q/^q = 
^fluid^^cxact _ I each combination of m/q and k/C, 
we show both '^Q^'"'\k,q,x) (top), and |A^'o/^'o| (bot- 
tom) with /max = 1 (solid lines), 2 (dotted lines) and 3 
(dashed lines). 

At large-scale {k/C = 10 or /c ~ 0.1 /i Mpc~^), neutri- 
nos with m/q = 10 become non-relativistic at around 
the horizon re-entry, and neutrinos with m/q = 100 
become non-relativistic well before the horizon re-entry 
{xnr = kTnr = 0.98 and 0.10 for m/q = 10 and 100, 
respectively), and fractional errors of the fluid approxi- 
mation peak at X = xh ~ 1. At small-scale {k/C = 100 
or fc ~ 1.0 h Mpc^^), neutrinos with m/q = 100 become 
non-relativistic at around the horizon re-entry, while neu- 
trinos with m/q = 10 become non-relativistic well af- 
ter the horizon re-entry {xnr = kr^r = 8.3 and 0.98 for 
m/q = 10 and 100, respectively). For neutrinos becom- 
ing non-relativistic after the horizon re-entry, fractional 
errors of the fluid approximation peak at x = Xnr- 



2. ■^i(k,q,x) 

Figure [3] shows the evolution of ^i{k, q, x) for two dif- 
ferent scales {k/C = 10 and 100) and two different mo- 
menta {m/q = 10 and 100). We calculate ^'f"''^(fc, g, x) 
and the fractional error in the fluid approximation in the 
same ways as before, while we calculate '^'^^'^^'^{k^q^x) 
from the exact solution of the Boltzmann equations given 
in Ea. (|^51) . Results are almost the same as the case 
for \E'o(A:, x), except for the asymptotic growth rate 
is ^"^"^^^{k^q^x) oc x^ cx a^/^ (see Eq.fgTI)). Since 
b{k, x) oc \E'o(fc, 9, a;) (X x^ and Q{k, x) oc ^\{k, g, x) oc x'^, 
the late-time evolution of bik^x) and 9{k,x) from the 
Boltzmann equations is consistent with the continuity 
equation (i.e., S{k,T) = —9{k,T)). 

3. ^2{k,q,x) 

Figure m shows the ratio of '^2{k, q, x) to '^o{k, q, x) for 
exact solutions (solid lines), /max = 2 (dotted lines) and 
^max = 3 (dashed lines) for two different scales {k/C = 10 
and 100) and two different momenta {m/q = 10 and 
100). As we discussed in § IIIIl we have a useful rela- 
tion between ^Q{k,q,x) and '^2{k, q, x) if we truncate 
the Boltzmann equations at /max = 2 (see Eq. ([M)l ). Here, 
we check how late-time evolution can simplify the rela- 
tion between 4'o(/e, g,x) and \l/2(fc, (?, x). At sufficiently 



^ The original paper did not take this factor of a^q/ao into 
account, and therefore, neutrino mass was overestimated by 
o.o/o,eq = 3000. This correction does not affect any qualita- 
tive/quantitative argument for ^i, but changes the interpreta- 
tions of the resulting and 9^ . 



Even though we do not include CDM explicitly, by setting 
(^ = '0 = 0, we are including CDM as a dominant source of 
the gravitational potential in the Boltzmann equations (i.e., the 
last term of Ea. (ll9l l). If we calculate tp{k,x) including the sup- 
pression of the gravitational potential due to massive neutrino 
free-streaming, the neutrino catch-up will be slightly slower. 
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^° 1 1 1 00 1 000 1 10 1 00 1 000 

X = kT X = kT 



FIG. 2: We show ^o{k, q, x) as functions oi x = kr with two different scales (fc/C = 10 and 100, where C = 0.0092 h Mpc~^) 
and two different momenta {m/q = 10 and 100). \E'o(fc,g, a;) is calculated from the exact solution, and the fractional difference 
is given as A^'o/'I'o = ^'q""*/5'q''^'^* — 1, where the solid line is for Zmax ~ 1, the dotted line is for /max ~ 2 and the dashed line 

is for Zmax = 3. 



late time, we have an asymptotic value of 'i!2{k^ q, x) = 
— |\['o(fc, q, x), which allows us to replace the anisotropic 
stress term in the fluid equations by the sound speed 
fEg . ([57]) l . As we see in the figure, the asymptotic value 
oi ^2{k,q,x)/^o{k,q,x) = —2/5 is reached at relatively 
early times. 



4- Summary 

We found that the applicability of the fluid approxima- 
tion on the perturbed distribution function, \E'/(fc, g, x), 
depends crucially on our choice of wavenumber, fc, mo- 
mentum, q and mass, to, of particles. Generally speak- 
ing, wavenumber, fc, sets the time of the horizon cross- 
ing, an = cl(x ~ 1), and the momentum and mass of 
the particles give the epoch of the relativistic to non- 
relativistic transition, Onr q/m. As long as a given 
particle becomes non-relativistic before the horizon cross- 
ing, flnr < gh, the fluid approximation and the exact 
solution of ^i{k,q,x) agree to better than 1% accuracy. 



Re- writing Om- < gh in terms of m/q and fc/C, we have, 

(57) 



q_ C 1/C 

771 k 4:\k 



This condition can be easily satisfied for a large-scale 
mode, where C/k > 1 (/c < 0.01 h Mpc^^), and the 
condition can be satisfied in a small-scale mode, C//s -C 1 
(A;>0.01 /iMpc-i), if Y > ^■ 

So far, we discussed the validity of the fluid approxima- 
tion for 4';(fc,q', x) with a given momentum and a given 
wavenumber. However, quantities such as S(k,x) and 
0(k,x) are given as integrals of g,a;) over momen- 
tum space (0 < g < oo) weighted by fo{q) and appropri- 
ate powers of q. The unperturbed distribution function, 
fo{q), exponentially suppresses the population of nearly 
relativistic neutrinos with an exponent of — ^ (see 

Ea. (|48p ). For a given lower limit of neutrino mass, 
nil, > 0.05 eV, and the current temperature of neutrinos, 
Ty.o 1.9 K, we have > 340. Therefore, the popula- 
tion of relativistic neutrinos with ^ ^ 0-1 is negligible, 
and docs not affect calculations of S(k,x) or 9{k,x). 
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FIG. 3: Same as Figure[2]for ^i{k,q,x). 



B. Limitation of Fluid Approximation on 5,j{k,x) 

Here, we study the limitation of the fluid approxima- 
tion on the neutrino density contrast, S^{k,x), and the 
velocity divergence, 9u{k, x), with two different masses of 
neutrinos, m,^ = 0.05 and 0.5 eV. As we have seen in the 
previous section, for given k and ao/aeq, the ratio of the 
neutrino mass to the current temperature of neutrinos 
determines the validity of the fluid approximation. 

With the current temperature of neutrinos, Tj^ q ^ 
1.9 K, we have m^/Tv,cq = 305 and 3050 for rrii, = 0.05 
and 0.5 eV, respectively. To find 5^{k,x) and 9i,(k,x), 
we integrate ^i{k,q,x) using Eas. (|2T|) ^ ((26l) for several 

^max- 



1. S,,{k,x) and9y{k,x) with small mu 

Figure [5] shows the evolution of the fractional er- 
rors of the fiuid approximation, /S.5u{k^x) / 5p{k^x) = 
5^'^'"^{k,x)/5l'''^''\k,x) ~ 1 and M^{k,x)/e^{k,x) = 
9^'^^'^{k,x)/0'^'^'^*'{k,x) — 1, as functions of a/ao for three 
different scales fc/C = 1, 10 and 100 with rrii, = 0.05 eV 
(C = 0.01 h Mpc-i). As expected, the fiuid approxi- 
mation does not yield accurate results for such a small 
mass. As was the case for ""^(fc, g, cc), the fractional 



error increases shortly after the horizon entrance, and 
then decreases as neutrinos become non-relativistic. For 
neutrinos with — 0.05 eV, Eq.© gives am/ao = 0.01, 
or Xnr = 9-4^, and as a result the fiuid approximation 
breaks down during entire evolution history of S^{k,x) 
and 9^{k, x) between the horizon crossing and the present 
epoch {uh < a < oq). Nevertheless, for the largest scale 
{k/C = 1), the error is below 10% level at low redshit, 
as the neutrinos become sufficiently non-relativistic. We 
also show the fractional errors of fluid approximation 
for (5^(fc,x) and 9^(k,x) calculated using the late time 
asymptotic value of Eq. dM]) : ^2{k,q,x) = —^^o{k,q,x) 
(dashed lines) . As we have seen in Figure IH this simple 
ansatz works well and follows the fractional error with 
^max = 2 at late time, a <C Onr- 



2. Sv{k,x) and8^{k,x) with large 

Figure [S] shows the evolution of the fractional errors 
of a fluid approximation as functions of a/ao for three 
different scales k/C = 1, 10 and 100 with rrii, = 0.5 eV. 
For neutrinos with rrii, = 0.5 eV, Eq.^ gives Onr/ao = 
0.001, or Xnr = 2.1^, and aU the scale with k/C > 0.5 
enter the horizon when neutrinos are relativistic. As a 
result, the fluid approximation with Zmax = 1 is only 
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1 100 10000 1 100 10000 

X = kT X = kT 



FIG. 4: We show 'ii2{k, q, x) /'iio{k, q, x) as functions of x = kr with two difTerent scales (k/C — 10 and 100) and two diflerent 
momenta {m/q = 10 and 100). Both ^o{k, q, x) and '^2(k,q,x) are calculated from the exact solution (solid line), or fluid 
approximation with /max = 2 (dotted line) and 3 (dashed line). 



accurate to ^ 1% at large scale {k/C ~ 1) and ^ 20% at 
small scale {k/C — 100) at low redshift. 

We see that Zmax = 2 and the ansatz approximates 
the small-scale density contrast and velocity divergence 
better than the case with lm_ax = 1, and the fluid ap- 
proximation becomes accurate to ^ 1% at large scale 
{k/C=l). 



C. Range of Validity of Fluid Approximation 

We have seen that the fractional errors of fluid approx- 
imation for (S,y(fc, x) and 9^{k, x) decrease for heavier par- 
ticles, but the errors are still significant on small-scales. 
Now, the question is what is the maximum wavenumber, 
kmax, below which we can use the fluid approximation 
with 10 or 20% accuracy for a given mass of neutrinos 
at a given time. Figure [7] shows the fractional error, 
/^5^{k,x)/5^{k,x) = 5l'^^'^{k,x)/5''J''^''\k,x) - 1, for four 
different masses of neutrinos at three different redshifts, 
= 0, 5 and 10. We find that the fluid approximation is 
only accurate to few~ 25% over a wide range of k at low 
redshift. 

Table [T] shows kmax with Zmax = 1 and 2 for various 



neutrino masses, — 0.05, 0.10, 0.50 and 1.0 eV, at 
flvc redshifts, z = 0, 1, 3, 5 and 10. The smaller the 
redshift is and the larger rrii, is, the larger kmax becomes. 
We sec that kmax is 3 ~ 4 times larger with Zmax = 2 
than with /max = 1- 

We are particularly interested in the fc-range of linear 
to mildly non-linear regime on the matter density power 
spectrum at the redshifts relevant to the future and on- 
going galaxy redshift surveys (z < 3), and for that pur- 
pose, 0.1 ;< kmax ^ 0.4 /i Mpc~^ will be necessary with 
sufficient accuracy: on smaller scales, the non-linearity is 
too large for power spectrum to be used for cosmology. 



VI. DISCUSSIONS AND CONCLUSIONS 

We have calculated the evolution of the perturbed 
distribution fimctions of massive neutrinos, 5';(fc, g, x), 
using the fluid approximation, i.e., truncation of the 
Boltzmann equations at /max = 1, 2 and 3. We com- 
pared the approximate solutions to the exact solutions 
that we have derived in this paper. When the distribu- 
tion function is dominated by the relativistic neutrinos, 
fluid approximation poorly represents the exact oscilla- 
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FIG. 5: (left): Time evolution of the fractional errors of A5^(fe, a;)/(5^(fc, a;) = 5^''''^{k, x) /ST'^'Hk, x) - 1. The solid lines 
show iniax = 1, while the dotted lines show imax ~ 2. The dashed lines show imax ~ 1, but with the ansatz for I — 2, 
^2{k,q, x) = -|*o(fc,g,2;). {right): Time evolution of the fractional errors of Ae,y{k,x)/6v(k, x) = 9^'^"^{k,x)/6''J''"'^ {k, x) - 1. 
Here, we use — 0.05 eV, and show the results at three different scales, k/C = 1, 10 and 100, corresponding to fc ~ 0.01, 0.1 
and 1.0 h Mpc~^, respectively. The vertical lines show the time of the horizon crossing for each mode. The present-day scale 
factor is oo = 3000. 



tion phase of '^i{k,q,x) calculated from the exact solu- 
tion. When the distribution function is dominated by 
non-relativistic neutrinos, \E'o(A:, g, x) and '^i{k,q,x) are 
sourccd mainly by the gravitational potential, ^{k,x), 
and decoupled from the higher multipoles, 5';>2(fc, x). 
This allovirs the fluid approximation to be an excellent ap- 
proximation to the growth of ^'o(A:,g,x) and 5'i(fc,q, a;) 
for small q. Then, we integrated the perturbed distribu- 
tion functions to calculate the quantities such as 8^(k, x) 
and 6u{k,x). Comparing the density contrasts of mas- 
sive neutrinos calculated from the fluid approximation 
to the exact solutions, we found that the fluid approxi- 
mation is only accurate to few^ 25% for k < 0.4 h Mpc~^ 
and 0.05 < m,j < 0.5eV To increase the accuracy of the 
fluid approximation further, it is necessary to either di- 
rectly solve for the Boltzmann hierarchy with Imax > 3 as 
in Eqs.([lH])--(lini), or solve fluid equations, Eqs.dH]) and 



(|28p. with an ansatz for an anisotropic stress, k^a{k,T), 
as we did for Imax = 2 in Eg. ([57]). 

We solved the Boltzmann equation for massive neutri- 
nos in an EdS universe for which (f){k,x) = ijj{k,x) = 0. 
In a more realistic multi-component fluid case, we have 
(j){k,x) ^ and ip{k,x) ^ due to the effect of mas- 
sive neutrinos even during the matter dominated epoch. 
Including this effect is straightforward. For the obser- 
vationally allowed range of neutrino masses, we expect 
the correction to be small, as the dominant source of 
gravitational potential is still CDM {f^, < 0.05). There- 
fore, our conclusions regarding the limitation of the fluid 
approximation is not affected by our using an EdS uni- 
verse. During the dark energy dominated epoch, the EdS 
approximation breaks down, and (j) ^-nd ip evolve. Nev- 
ertheless, the correction will be limited to the dark en- 
ergy dominated epoch, a > ajjE, and the scale around 
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FIG. 6: Same as Figure [5] for = 0.5 eV. 



^max — 1 


= 0.05 eV 


0.1 eV 


0.5 eV 


1.0 eV 


z=0 
z=l 
z=3 
z=5 
z=10 


0.009 (0.032) 
0.007 (0.021) 
0.005 (0.015) 
0.003 (0.012) 
0.002 (0.009) 


0.018 (0.064) 
0.012 (0.042) 
0.009 (0.028) 
0.008 (0.023) 
0.004 (0.017) 


0.090 (0.31) 
0.060 (0.19) 
0.039 (0.13) 
0.032 (0.097) 
0.023 (0.068) 


0.18 (0.62) 
0.12 (0.38) 
0.079 (0.25) 
0.060 (0.19) 
0.042 (0.14) 


^max — 2 


= 0.05 eV 


0.1 eV 


0.5 eV 


1.0 eV 


z=0 
z=l 
z=3 
z=5 
z=10 


0.037 (0.097) 
0.024 (0.064) 
0.016 (0.042) 
0.013 (0.034) 
0.009 (0.024) 


0.073 (0.19) 
0.045 (0.13) 
0.032 (0.079) 
0.024 (0.064) 
0.018 (0.045) 


0.36 (0.94) 
0.22 (0.58) 
0.15 (0.36) 
0.11 (0.29) 
0.079 (0.19) 


0.72 (1.88) 
0.44 (1.16) 
0.27 (0.72) 
0.22 (0.54) 
0.15 (0.38) 



TABLE I: The maximum wavenumber, k„iax[h Mpc for which the fluid approximation is accurate at 10 (20)% or better. 



k ^ kYs,{o.DE) ^ ^Fs(anr). We found that, as long 
as the term proportional to the gravitational potential, 
ip{k,x), dominates the right hand side of Ea. ([T5)) . the 
fluid approximation is valid. Therefore, unless the effect 



of the dark energy suppresses ip^k^x) much faster than 
the growth of e^{q,x) oc a^, k^ax at a > ude will not 
change significantly. 

Since we have studied the evolution of the distribution 
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FIG. 7: We show the fractional errors, M^{k,x)/&u{k,x) = 5^'^''^{k,x)/5T'"'\k,x) - 1, for four different masses of neutrino at 
three different redshifts, z = 0, 5 and 10 as functions of wavenumber, where the thick and thin hnes are for Zmax = 1 and 2, 
respectively. 



function solving the collision-less Boltzmann equation, 
one can apply these results to other collision-less particles 
in general. 

Now, why is fluid approximation useful? Future and 
on-going dark energy missions aim at the accurate mea- 
surement of the galaxy /matter power spectrum with an 
accuracy better than 1%. One might think that the cos- 
mological linear perturbation theory has already been 
well established, and the numerical codes such as CMB- 
fast and CAMB can calculate the linear matter power 
spectrum with an accuracy better than 1%. 

However, the linear perturbation theory breaks down 
at small-scale and low redshift, where the density con- 
trast becomes non-linear (fc > 0.1 h Mpc^^ at z ~ 1) 
[314^. Therefore, in order to exploit the cosmological 
information contained in a given survey, one needs to un- 
derstand the non-linearities on the galaxy/matter power 
spectrum [sol - fs^ . 

Among the non-linearities, the matter clustering has 
been well understood in the mildly non-linear regime (see 
[ssj , for a review), but the theories have been limited to 
CDM dominated universe. The pressure gradient term 
in the Euler equation was completely ignored. 

In our previous work, we developed the 3rd-order per- 



turbation theory with the pressure gradient terms explic- 
itly included [ij] (also see [53 - l57| ). With this extension 
to the higher order perturbation theory as well as within 
the limitation on the accuracy of (5^ calculated from the 
fluid approximation, we can now calculate the next-to- 
linear order matter power spectrum with massive neu- 
trino free-streaming effect, properly included. 

Since the structure formation is mostly affected by the 
most massive species of neutrinos, and the current con- 
straints on the total mass of neutrinos indicate that at 
least one of the neutrino species has a mass of order a 
tenth of eV, the use of fluid approximation is limited with 
an accuracy of few to 25% over k < 0.4 h Mpc~^ for 
z < 10. As a result, for a small fraction of massive neu- 
trino, fi, < 0.04 for ^ • m^^i < 0.5 eV, the fractional error 
on the matter density contrast, (5m = (1 — fv)5c + fu^v^ 
calculated with the fluid approximation is accurate to 
sub-percent level. 

This material is based in part upon work supported by 
the Texas Advanced Research Program under Grant No. 
003658-0005-2006, by NASA grants NNX08AM29G and 
NNX08AL43G, and by NSF grant AST-0807649. M. S. 
thanks for warm hospitality of Astronomical Institute at 
Tohoku University where part of this work was done. 
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Appendix A: Sound speed versus Velocity Dispersion 

The wavenumber corresponding to the free-streaming scale, /cpS; is defined by the single- fluid continuity and Euler 
equations: 

(j(k,T)+0(k,T) =0 (Al) 

^(k, r) + H{T)9{k, r) + cI{t) [klsir) - k^r)] 5(k, r) = 0, 

(A2) 

where 

fcps(r) = ^/|^, (A3) 
V 2 Cs(t) 

is the scale that divides the characteristics of the solution, (5(k, r). For k < kps, ^(k, r) grows, while for k > fcps, 
(5(k, r) oscillates. 

In the literature, however, the sound speed, Cs(t), has often been replaced, or crudely approximated, by the velocity 
dispersion, CTi/, such that Cg — Uiyj without any justification. Strictly speaking, the velocity dispersion should not be 
used to define the free-streaming scale, /cpSj as the Euler equation contains sound horizon, Cg = not the velocity 
dispersion. In this Appendix, we shall clarify this issue. 

The sound speed is given by 

. SP{k,r) _ 1 Iq'dqj£^fo{q)^o{k,q,r) 
' > - Sp{k,T) 3jq^dqeiq,r)Mq)^o{k,q,Ty ^ ' 

and the velocity dispersion is given by 



Sq^dq 


q 




J q^dqfoiq) 



J q^dqfo{q) 

Note that the sound speed can depend on k in general, while velocity dispersion is, by definition, independent of k. 
As we see from Eq. (|46l) . in the non-relativistic limit, where the phase-space distribution of the particles is dominated 



by non-relativistic particles {q <^ a{x)m), the k dependence of 4'o(fc,q, r) = ^o(fc, ''") in'g'^'' ^^^^ canceled both in 
the denominator and numerator of Ea. (IA4l) . If perturbations arc adiabatic, i.e., ^ = "^^^7' have 

, , , P(t) , , wfr) , . 

p(r) Zn{T){\ ^ w(t)\ 

where w{t^ = is an equation of state. Since the velocity dispersion in the non-relativistic limit is given as 



2/ , , 1 J q^dq q^fo{q)'io{k,q,T) 
a^{T)m? j q^dqJo{q)-^o{k,q,T) 



we have 

"''^'^ 3"''^'^^9l + ia3(r) - 9 



cl{r) ^ + l-^^4rT - l-lir)- (A8) 



Here, we have used aKr) <C 1. Therefore, in the non-relativistic limit, we have = — 0.745cry,i. 

Appendix B: Exact Analytic Solution of '^i(k,q,x) 
We derive the analytic solutions of the Boltzmann equation. We start from Eq.([T 



d'^ik,q,fi,x) . q , . , . 
T, hi-^ -fi-9{k,q,fj,,x) = b{k,q,fj,,x), (Bl) 

ox £[q,x) 



where 



c,/, \ - ■'^ili^) in \ 9(t){k,x) 

S(k,q,^i,x) = I ^iipik.x) . 

q ax 

The solution of the first order linear differential equation with the source term S{k^ q, /i, x) is 



^{k, q, fj., x) = q, 11, Xi) cxp 



X / dx' cxp ifi / 

J Xi J Xi 



e{q,x") 



dx" 



eiq,x) 
S{k, q, ^, x') 



-dx' 



cxp 



-dx' 



where we define 



z{x) 



c 



e{q,x 



-dx', 



for arbitrary constant, C E TZ. 
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(B2) 



(B3) 



(B4) 



1. massless case with constant ((> and ^ 

For massless neutrinos, e(q, x) = q, with constant gravitational potential, (j) = ~ 0, we can avoid the complexity 
in the time dependent source term, S{q, x), and the form of Eq. (jB3p is simplified to 

^{k,q,fi,x) = §(fc,g,/i,xOe-'^("^-"^') / dx'e-*^(="-^ ) 



^ — ifj.{x~Xi) 



^(fc). 



(B5) 



= '^{k■,q,^i■,x^) -ipik) 

Now, we expand Eq. ljBSp by Legcndrc polynomials, using the following identities and orthogonality condition, 

q, 11, x) = ^(-*)'(2/ + q, x)Pi{p), (B6) 

e-^>^^-~^"^ = Y.{-i)\2l + l)3i{x - x')Pi{fi), (B7) 

2 



j ^diJiPi{ti)Pi.{ti) 



21 + 1 



-5w. 



(B8) 



We find 



Y,{-iY{2l' + l)^l,{k,q,x)Pv{^i) = ^^(-z)''+'"(2/' + l)(2r + l)^i,{k,q,x,)jv.{x - x,)P,{^l)Pri^l) 



I' I" 



(B9) 



Multiplying the both sides by and integrating over /x, we find. 



/ I' I" 




3i0 



(BIO) 
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where we have used the Wigner 3-j symbols to write 

"1 dn 



(BU) 



Therefore, for a given I, we have 



*,(fc,9,a:)=^^(-i)''+'''-'(2r + l)(2r + l)*i,(fc,g,x,)>'(a:-x,) (^J^ ) +iPik) [Sm - jiix ^ x,)] .{BU) 



Here, non-zero Wigner 3-j symbols must satisfy the triangular inequalities such that 

\l-l'\ < I" <l + l', 

and from the initial conditions, we have ^^I'ysixi) ~ 0. 
The exact solutions for Z = and 1 are 



(B13) 



^o{k,q,x) = ^o{k,q,Xi)jo{x - Xi) -3'^i{k,q,Xi)ji{x - Xi) 

+ 5^2ik,q,x,)j2{x - X,) +^{k)[l - ja{x - X,)], (B14) 

^i{k,q,x) = ^o{k,q,Xi)ji{x - Xi) + ^i{k,q,Xi)ja{x - Xi) - 2'^2ik,q,Xi)ji{x - Xi) 

- 2^i{k,q,Xi)j2{x - Xi) + 2,^2(k,q,Xi)jz{x - Xi) - ip{k)ji{x - Xi). (B15) 

2. massive case with constant and ?/> 

For massive neutrinos, e{q,x) ^ q, with constant gravitational potential, = -0 = 0, we have 



(B16) 



where we have used Eqs. (|B6|) ^ (jB8p and z — z' = z{x) — z{x') = J^, dx' j^-^^ as defined in Ea. (jB4 
Multiplying the both sides by Pi{ij) and integrating over /i, we find 



Y,{-iy^l{k,q,x) EE(-*)''^'"(2^' + l)(2^" + l)*;'(fc,'Z,a;,)ji"(^-^0 
I I' I" 



I V l"\ 



(B17) 



Therefore, for a given Z, we have 



g, x) = Y.Y.^-iY+'"-\2l' + 1)(2Z" + 1)^-^ (fc, g, - z,) 



V I" 



I V I" 




ip{k) / dx 



,e{q,x' 



2/ + 1 



2Z + 1- 



(B18) 



The exact solutions for / = and 1 are 



'^o{k,q,x) = ^Q{k,q,Xi)ja{z - z.,) - 3*i(fc, g, (z - z,) + 5'^2ik,q,Xi)j2iz - z,) 

+ i^ik) r dx'^-^^n{z'z'), 

^i(k,q,x) = ^!o{k,q,Xi)ji{z - Zj) + *i (fc, g, Xj)jo(z; - z^) ~ 2'^2{k,q,Xi)ji(z - Zi) 

- 2^i(fc, g,.T,;)j2(z - Zi) + 'i^2{k,q,Xi)jz{z - z^) 



0(fc) r dx'^^^^ 



1 2 

-^joiz ~ z') - -j2{z - z') 



(B19) 



(B20) 
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3. general case 



For more general cases, where neutrinos are either massless or massive, and 0^0 and ip ^ Q. In this case, we use 



the full expression of the source term, S{k, q, fi, x) = i^-^^^^^ip{k, x) 



^,^j3fl^^^k,x)-^^. 

Again, we expand Ea. (jB3[) with the scries of Lcgendre polynomials with the time dependent source term, and find 



I' 



+ iM5](-*)''(2/' + l)Pr(M) r dx'^-^^{k,x')jv{z-z') 



^(-^)''(2/' + 1)P.(m) fdx'^-^^Hz - z'). 



(B21) 



Multiplying the both sides by and integrating over /x, we find 



Y,{-i)'^i (fc, g, .t) = ^ ^(-*)''+'" {21' + 1)(2/" + I)*,, (fc, g, x^Ji" - z.,) 



V I" 



J2{^ty r dx''-^^ij{k,x') 

^ .w r , , d(l){k,x') . 



I 



I V I" 
0^ 

l + l 



21 + 1 



(B22) 



With the recursion relation. 



we have 



2^ + 1 



/ . , , l+l . . . 



2; + 1" 



(B23) 



^lik, q, x) = J2 Y.{-iY+'"~\2l' + l)(2r' + 1)'!',, (fc, q, - z,) 



V I" 



I V I"' 
,0 0, 



dx' 



eiq,x') 



^{k,x')~-^cjy{k,x') 



q eiq,x') 
+ (f>{k,Xi)ji{z - Zi) ~ (l){k,x)Sio, 



I . I ,^ ^ + 1 

Ji-i{z - Z ) 



2/ + 1 



21 + 1 



ji+i{z - z') 



(B24) 



for a fixed I. We can easily recover the massless and massive case solutions Eq. (jB12[ ) and (jBlSp from Eq. (jB24[ ) with 
approximations such as 'ijj{k,x) = (f>{k,x) = and/or e{q,x) = q. 
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